﻿module FreeEnergy
    // Still need to make the Poly1,2,3 functions more type safe
    let unknown = '?'
    let unknown2 = -1
    let IssueError msg = 
        printfn "%s" msg



    // Reads a triangle/symmetric matrix listing all the bond strengths between different glues
    let ConvertArray2DInt (a : int array) width length = 
        let M = Array2D.init width length (fun i j -> a.[i * length + j])
        M

    let ConvertArray2DChar (a : char array) width length = 
        let M = Array2D.init width length (fun i j -> a.[i * length + j])
        M

    let getFE (a : int[]) (glueMatrix : int[,]) =
        let mutable s = 0
        for i = 0 to a.Length - 2 do
            s <- s + glueMatrix.[a.[i], a.[i+1]]
        s

    let getMaxLocalFE (a : int[,]) (glueMatrix : int[,]) =
        Array.init (a.GetLength(0)) (fun i -> Array.init (a.GetLength(1)) (fun j -> a.[i, j]))
            |> Seq.map (fun i -> getFE i glueMatrix) |> Seq.max

    let getMinLocalFE (a : int[,]) (glueMatrix : int[,]) =
        Array.init (a.GetLength(0)) (fun i -> Array.init (a.GetLength(1)) (fun j -> a.[i, j]))
            |> Seq.map (fun i -> getFE i glueMatrix) |> Seq.min

    let getMaxFE (glueMatrix : int[,]) = 
        Array.init (glueMatrix.GetLength(0)) (fun i -> Array.init (glueMatrix.GetLength(1)) (fun j -> glueMatrix.[i, j]))
            |> Seq.map (fun i -> i |> Seq.max) |> Seq.max

    let getMinFE (glueMatrix : int[,]) =
        (Array.init (glueMatrix.GetLength(0)) (fun i -> Array.init (glueMatrix.GetLength(1)) (fun j -> glueMatrix.[i, j]))
            |> Seq.map (fun i -> i |> Seq.min) |> Seq.min)


    // Get the difference between the max and min of all the entries in glueMatrix
    let getD (glueMatrix : int[,]) = 
        (getMaxFE glueMatrix) - (getMinFE glueMatrix)

    let convertToBase (x : int) (b :int) =
        let mutable temp = x
        //let mutable step = b
        let mutable index = 0
        let mutable A = [| |]
        while temp > 0 do
           A <- Array.append A [| temp % b |]
           index <- index + 1
           temp <- temp / b
        done
        A


    let CountFreePairs (a : int[]) (glueMatrix : int[,]) =
        let mutable sum = 0
        let b = glueMatrix.GetLength(0)
        let unknowns = a |> Seq.map (fun i -> if i = unknown2 then 1 else 0) |> Seq.sum
        
        // Set trialCount to b ** unknowns
        let mutable trialCount = 1
        for i = 1 to unknowns do
            trialCount <- trialCount * b
        done
        
        let genBinary = (System.Random()).Next(2)
        let minFE = getMinFE glueMatrix
        let mutable polynomial = Array.init ((getMaxFE glueMatrix) - minFE + 1) (fun i -> 0)

        // Goes through all the permutations that the string can take
        for i = 0 to trialCount - 1 do
            // Takes the trial and converts it into the array for the changes that will be applied
            let permutationTemp = convertToBase i b
            
            // Pad the array to the length specified
            let permutation = Array.append permutationTemp (Array.init (a.Length - permutationTemp.Length) (fun i -> 0))
            let mutable stackIndex = 0
            let mutable trial = [| |]

            // Then applies that permutation to the string
            for j = 0 to a.Length - 1 do
                if a.[j] = unknown2 then
                    trial <- Array.append trial [| permutation.[stackIndex] |]
                    stackIndex <- stackIndex + 1
                else
                    trial <- Array.append trial [| a.[j] |]

            polynomial.[(getFE trial glueMatrix) - minFE] <- polynomial.[(getFE trial glueMatrix) - minFE] + 1
        polynomial

    // Gives an eleme
    let Poly1 (el : int) (a : int) (b : int) (glueMatrix : int[,]) =
        if el < 1 then
            IssueError "The length of the poly is less than 1"
            [| |]
        else
            let template : int[] =
                if el = 1 then
                    [| a |]
                else
                    Array.append [| a |] (Array.init (el - 2) (fun i -> unknown2)) |> Array.append <| [| b |]
            
            CountFreePairs template glueMatrix

    let Poly2 (el : int) (a : int) (b : int) (d1 : int) (d2 : int) (glueMatrix : int[,]) =
        if el < 1 then
            IssueError "The length of the poly is less than 1"
            [| |]
        else
            let template : int[] =
                if el = 1 then
                    [| a |]
                else
                    Array.append
                        (Array.append [| a |] (Array.init (el / 2 - 2) (fun i -> unknown2)) |> Array.append <| [| d1 |])
                        (Array.append [| d2 |] (Array.init (el / 2 - 2) (fun i -> unknown2)) |> Array.append <| [| b |])
            
            CountFreePairs template glueMatrix



    let Poly3 (el : int) (a : int) (b : int) (d1 : int) (d3 : int) (d2 : int) (glueMatrix : int[,]) =
        if el < 1 then
            IssueError "The length of the poly is less than 1"
            [| |]
        else
            let template : int[] =
                if el = 1 then
                    [| a |]
                else
                    Array.append
                        (Array.append (Array.append [| a |] (Array.init (el / 2 - 2) (fun i -> unknown2)) |> Array.append <| [| d1 |])
                            [| d3 |])
                        (Array.append [| d2 |] (Array.init (el / 2 - 2) (fun i -> unknown2)) |> Array.append <| [| b |])
            
            CountFreePairs template glueMatrix

    let PolyMultiply (a : int[]) (b : int[]) =
        let mutable result = Array.init (a.Length + b.Length - 1) (fun i -> 0)
        for i = 0 to a.Length - 1 do
            for j = 0 to b.Length - 1 do
                result.[i + j] <- result.[i + j] + a.[i] * b.[j]
        result

(*    let rec Build el = 
        if el = 1 then
        elif el > 1 then
            let mutable list = Build (el / 2) // el divided by 2 is floored
            if el % 2 = 0
        else
            IssueError "el must be greater than or equal to 1"
            -1

    let rec Build (el : int) (glueSet : int array) = 
        if el < 1 then
            IssueError "el must be take on positive values"
        elif el = 1 then
            let mutable function3 el a b d1 d3 d2 = null
            let mutable function1 el a b = if a = b then 1 else 0
            [| (function1, function3) |]
        else
            let mutable LIST = Build (el / 2) glueset // el is floored
            if (el % 2) = 0 // el is even
                let x = 0
                let function1 
                let function2 el a b d1 d2 <- PolyMultiply   <| PolyMultiply (funtion1 (el / 2) a d1) (function1 (el / 2) d2)
            else // el is odd
        
            

    let rec BuildRecursive el = 
        if el > 2
            
        else
            
        done*)
            
    (*let RecursiveExtract = 
        [| |]

    let Extract el e (glueMatrix : int[,]) = 
        let mutable a = -1
        let mutable b = -1
        for i = 0 to glueMatrix.GetLength(0) - 1 do
            for j = 0 to glueMatrix.GetLength(1) - 1 do
                let poly = Poly1 el i j glueMatrix
                if poly.[e] <> 0 then
                    a <- i
                    b <- j
            done
        done
        if el = 1 then
            [| |]
        else
            RecursiveExtract*)


    let l_star (n : int) (k1 : int) (k4 : int) (delta : float) =
        let c1 = 2.0 + delta
        let c_star = c1 / ((c1 - 2.0) * log 2.0)
        let c2 = (c1 / 2.0) * ((log c_star) / (log 2.0)) + 2.5 - (1.0 / (log 2.0))
        ceil(c1 * ((log <| float n) / (log 2.0)) + c2 * (float <| max k1 k4))



    // Replaces all the entries in matrix a, with the index they have in glueSet
    let ConvertToIntMatrix (a : char[,]) (glueSet : char[]) =
        let convert = Array.init 256 (fun i ->
            let mutable x = -1
            for j = 0 to glueSet.Length - 1 do
                if i = int glueSet.[j] then
                    //printfn "%c;%d" (glueSet.[j]) i
                    x <- j
            x)
        Array2D.init (a.GetLength(0)) (a.GetLength(1)) (fun i j -> convert.[int a.[i, j]])



    let SpecialConcatenate (x : int[]) (y : int[]) = 
        Array.append y.[0 .. ((y.Length / 2) - 1)] x |> Array.append <| y.[(y.Length / 2) .. (y.Length - 1)]

    // Satisfies constraint C1, C2, C3, C4, C5, C6, C9 for DNA glueSets
    let FastDWD (delta : float) n k1 k2 k3 k4 k5 k6 (glueMatrix : int[,]) = 
        let glueSet = [| 'A'; 'C'; 'G'; 'T' |]
        if (glueMatrix.GetLength(0) <> 4) && (glueMatrix.GetLength(1) <> 4) then
            IssueError "'glueMatrix must be a 4 by 4 matrix for DNA bases"
            Array2D.init 0 0 (fun i j -> unknown)
        else
            let el = ExpCount.findOptimalL n (int <| l_star (int n) (int k1) (int k4) (float delta)) k1 k4
            let Mchar = dnaTools.FastDWD1to6 n el k1 k2 k3 k4 k5 k6
            let Mint = ConvertToIntMatrix Mchar glueSet
            let Wmax = getMaxLocalFE Mint glueMatrix
            let Wmin = getMinLocalFE Mint glueMatrix
            printfn "%d:%d=%d" Wmax Wmin (3* (getD glueMatrix))

            if (Wmax - Wmin) <= 3 * (getD glueMatrix) then
                Mchar
            else
                let L = el * 2
                let b = glueMatrix.GetLength(0)

                // Set trialCount to b ** unknowns
                let mutable trialCount = 1
                for i = 1 to L do
                    trialCount <- trialCount * b
                done
        
                let minFE = getMinFE glueMatrix
                let mutable polynomial = Array.init ((getMaxFE glueMatrix) - minFE + 1) (fun i -> 0)
        
                let W = Array.init (Mint.GetLength(0)) (fun i -> Array.init (Mint.GetLength(1)) (fun j -> Mint.[i, j]))
                let mutable W2 = Array2D.init (Mint.GetLength(0)) L (fun i j -> '?')
                
                let pad (a : int[]) (desiredLength : int) =
                    Array.append a (Array.init (desiredLength - a.Length) (fun i -> 0))


                let mutable previous = pad (convertToBase 0 b) L
                let mutable current = pad (convertToBase 1 b) L
                let mutable freeEnergies = [| (getFE current glueMatrix); (getFE previous glueMatrix) |]
                let mutable delta = abs <| (getFE current glueMatrix) - (getFE previous glueMatrix)
                
                // Goes through all the permutations that the string can take
                for i = 2 to trialCount - 1 do
                    previous <- current
                    current <- pad (convertToBase i b) L
                    freeEnergies <- Array.append freeEnergies [| getFE current glueMatrix |]
                    delta <- max delta <| (getFE current glueMatrix) - (getFE previous glueMatrix)
                done
        
                for i = 0 to W.Length - 1 do
                    let mutable toAdd = [| |]
                    let mutable j = trialCount - 1
                    let mutable cond = true
                    while (j >= 0) && cond do
                        j <- j - 1
                        if (abs (freeEnergies.[j] - (getFE W.[i] glueMatrix))) <= delta then
                            cond <- false
                            toAdd <- pad (convertToBase j b) L
                            printfn "WOOT"
                    done

                    if toAdd.Length > 0 then
                        let toAdd2 = SpecialConcatenate W.[i] toAdd
                        for j = 0 to L - 1 do
                            W2.[i, j] <- glueSet.[toAdd2.[j]]
                        done
                done
                printfn "%A\n\n\n" W
                W2